Introduction


Host plant genotype and phenotype influence arthropod community composition. They influence fitness, can mediate competition, predation and host-parasite interactions. Plant size, leaf chemistry, defense compounds, trichome density etc all affect species in different ways. However, it is hard to pin down general principles in how host plant characteristics affect herbivores across species Insect traits can be used to provide a mechanistic connection between species and environment. Insect traits have increasingly been used to understand insect response to characteristics of landscape and understand community composition.


I set out with the following objectives:

(i) How do herbivore communities respond to phenotypic variation?
(ii) Which plant traits account for herbivore community responses to host-plant phenotype?
(iii) Which insect traits account for herbivore community responses to host-plant phenotype




Results



Some context: Richness across hybrid sites



(i) How do herbivore communities respond to phenotypic variation?


p<-metaMDS(hem.both, distance="bray", k=2, trymax=100)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2614694 
## Run 1 stress 0.2591777 
## ... New best solution
## ... Procrustes: rmse 0.07578497  max resid 0.2182768 
## Run 2 stress 0.2643395 
## Run 3 stress 0.2711564 
## Run 4 stress 0.2591012 
## ... New best solution
## ... Procrustes: rmse 0.004049879  max resid 0.02030985 
## Run 5 stress 0.2680823 
## Run 6 stress 0.2669241 
## Run 7 stress 0.2676229 
## Run 8 stress 0.2668678 
## Run 9 stress 0.2686223 
## Run 10 stress 0.264573 
## Run 11 stress 0.2650065 
## Run 12 stress 0.2764339 
## Run 13 stress 0.3046596 
## Run 14 stress 0.2603386 
## Run 15 stress 0.2646121 
## Run 16 stress 0.2642043 
## Run 17 stress 0.260578 
## Run 18 stress 0.264796 
## Run 19 stress 0.2646198 
## Run 20 stress 0.2642683 
## Run 21 stress 0.2708115 
## Run 22 stress 0.261162 
## Run 23 stress 0.2654999 
## Run 24 stress 0.2630344 
## Run 25 stress 0.2657158 
## Run 26 stress 0.2587648 
## ... New best solution
## ... Procrustes: rmse 0.01136149  max resid 0.04027083 
## Run 27 stress 0.2617533 
## Run 28 stress 0.263632 
## Run 29 stress 0.2749018 
## Run 30 stress 0.2654376 
## Run 31 stress 0.2621657 
## Run 32 stress 0.2667588 
## Run 33 stress 0.265062 
## Run 34 stress 0.2766734 
## Run 35 stress 0.2720742 
## Run 36 stress 0.2680699 
## Run 37 stress 0.2671735 
## Run 38 stress 0.2703416 
## Run 39 stress 0.2651467 
## Run 40 stress 0.2658195 
## Run 41 stress 0.2629638 
## Run 42 stress 0.2648128 
## Run 43 stress 0.2630296 
## Run 44 stress 0.2613272 
## Run 45 stress 0.2648603 
## Run 46 stress 0.2618145 
## Run 47 stress 0.27441 
## Run 48 stress 0.267981 
## Run 49 stress 0.266718 
## Run 50 stress 0.2606572 
## Run 51 stress 0.2696575 
## Run 52 stress 0.2667685 
## Run 53 stress 0.2636413 
## Run 54 stress 0.2686025 
## Run 55 stress 0.2697989 
## Run 56 stress 0.2589156 
## ... Procrustes: rmse 0.009479284  max resid 0.04055719 
## Run 57 stress 0.2648442 
## Run 58 stress 0.2603029 
## Run 59 stress 0.2696071 
## Run 60 stress 0.2615379 
## Run 61 stress 0.2652593 
## Run 62 stress 0.2659076 
## Run 63 stress 0.2666821 
## Run 64 stress 0.2682568 
## Run 65 stress 0.2668895 
## Run 66 stress 0.2711279 
## Run 67 stress 0.2657382 
## Run 68 stress 0.2654936 
## Run 69 stress 0.2676306 
## Run 70 stress 0.2818346 
## Run 71 stress 0.2822666 
## Run 72 stress 0.2655272 
## Run 73 stress 0.2660869 
## Run 74 stress 0.264887 
## Run 75 stress 0.2619298 
## Run 76 stress 0.2634607 
## Run 77 stress 0.2650585 
## Run 78 stress 0.2708242 
## Run 79 stress 0.2644942 
## Run 80 stress 0.2731182 
## Run 81 stress 0.2611793 
## Run 82 stress 0.2629705 
## Run 83 stress 0.2653446 
## Run 84 stress 0.2772961 
## Run 85 stress 0.2698972 
## Run 86 stress 0.2651032 
## Run 87 stress 0.2745211 
## Run 88 stress 0.2655207 
## Run 89 stress 0.2652121 
## Run 90 stress 0.263763 
## Run 91 stress 0.2619298 
## Run 92 stress 0.2654345 
## Run 93 stress 0.2731397 
## Run 94 stress 0.2630328 
## Run 95 stress 0.2627176 
## Run 96 stress 0.2668713 
## Run 97 stress 0.2614056 
## Run 98 stress 0.2603036 
## Run 99 stress 0.265145 
## Run 100 stress 0.2633985 
## *** No convergence -- monoMDS stopping criteria:
##    100: stress ratio > sratmax
ordiplot(p,type="n")
orditorp(p,display="species",col="red",air=0.01)
orditorp(p,display="sites",cex=1.25,air=0.01)

ordiplot(p,type="n")
orditorp(p,display="species",col="red",air=0.01)
orditorp(p,display="sites",pch=19, col="blue", label=F)
ordiellipse(p,groups=env.both$morphotype,draw="polygon",col="grey90",label=T, kind="se")



Are communities distinctly different between different morphotype trees within the same site?


hem.ero<-hem[env$site=="ERO",]
hem.ali<-hem[env$site=="ALI",]
hem.kaiy<-hem[env$site=="KAIY",]
hem.olaa<-hem[env$site=="OLAA",]
hem.kaio<-hem[env$site=="KAIO",]

#par(mfrow=c(2,2))

moment<-metaMDS(hem.ero, k=3)
## Wisconsin double standardization
## Run 0 stress 0.06827439 
## Run 1 stress 0.06827615 
## ... Procrustes: rmse 0.0004506555  max resid 0.0007439575 
## ... Similar to previous best
## Run 2 stress 0.06827984 
## ... Procrustes: rmse 0.0008459453  max resid 0.001391683 
## ... Similar to previous best
## Run 3 stress 0.068273 
## ... New best solution
## ... Procrustes: rmse 0.0003991526  max resid 0.0007228103 
## ... Similar to previous best
## Run 4 stress 0.08827491 
## Run 5 stress 0.0822194 
## Run 6 stress 0.08132868 
## Run 7 stress 0.08854191 
## Run 8 stress 0.06827149 
## ... New best solution
## ... Procrustes: rmse 0.003514993  max resid 0.006012255 
## ... Similar to previous best
## Run 9 stress 0.1413167 
## Run 10 stress 0.0882751 
## Run 11 stress 0.06827118 
## ... New best solution
## ... Procrustes: rmse 0.002863892  max resid 0.005034434 
## ... Similar to previous best
## Run 12 stress 0.08132942 
## Run 13 stress 0.06827011 
## ... New best solution
## ... Procrustes: rmse 0.0005564563  max resid 0.001017243 
## ... Similar to previous best
## Run 14 stress 0.08133052 
## Run 15 stress 0.08827486 
## Run 16 stress 0.2421808 
## Run 17 stress 0.06828125 
## ... Procrustes: rmse 0.00228561  max resid 0.003866028 
## ... Similar to previous best
## Run 18 stress 0.06827822 
## ... Procrustes: rmse 0.002354057  max resid 0.003952792 
## ... Similar to previous best
## Run 19 stress 0.06827094 
## ... Procrustes: rmse 0.002173734  max resid 0.003670135 
## ... Similar to previous best
## Run 20 stress 0.06827642 
## ... Procrustes: rmse 0.001993287  max resid 0.003361344 
## ... Similar to previous best
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main="Escape road old flow")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="ERO",]$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env[env$site=="ERO",]$morphotype == "P")

Lab<-c( "Hybrid", "Pubescent")
color<-c("blue", "red")
legend(  "topright", legend=Lab,pch=c(15,19),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

moment<-metaMDS(hem.kaiy, k=3)
## Wisconsin double standardization
## Run 0 stress 0.06676995 
## Run 1 stress 0.06676756 
## ... New best solution
## ... Procrustes: rmse 0.0004757995  max resid 0.000815727 
## ... Similar to previous best
## Run 2 stress 0.08024762 
## Run 3 stress 0.0802441 
## Run 4 stress 0.08845382 
## Run 5 stress 0.08024536 
## Run 6 stress 0.06676711 
## ... New best solution
## ... Procrustes: rmse 0.004773203  max resid 0.008547457 
## ... Similar to previous best
## Run 7 stress 0.08627771 
## Run 8 stress 0.08024788 
## Run 9 stress 0.08627932 
## Run 10 stress 0.0667668 
## ... New best solution
## ... Procrustes: rmse 0.0001037532  max resid 0.0002291509 
## ... Similar to previous best
## Run 11 stress 0.08848844 
## Run 12 stress 0.08847556 
## Run 13 stress 0.066772 
## ... Procrustes: rmse 0.0008988719  max resid 0.001495977 
## ... Similar to previous best
## Run 14 stress 0.08627253 
## Run 15 stress 0.06676287 
## ... New best solution
## ... Procrustes: rmse 0.001662798  max resid 0.002918259 
## ... Similar to previous best
## Run 16 stress 0.06676659 
## ... Procrustes: rmse 0.001537349  max resid 0.002704864 
## ... Similar to previous best
## Run 17 stress 0.06677093 
## ... Procrustes: rmse 0.003762415  max resid 0.006646404 
## ... Similar to previous best
## Run 18 stress 0.06676934 
## ... Procrustes: rmse 0.003467202  max resid 0.006118675 
## ... Similar to previous best
## Run 19 stress 0.06676371 
## ... Procrustes: rmse 0.0007865925  max resid 0.001268778 
## ... Similar to previous best
## Run 20 stress 0.06676702 
## ... Procrustes: rmse 0.002951224  max resid 0.005323313 
## ... Similar to previous best
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main="Kaiholena Young flow")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="KAIY",]$morphotype == "G")
points(mds.fig, "sites", pch = 19, col = "red", select = env[env$site=="KAIY",]$morphotype == "P")

Lab<-c("Glabrous",  "Pubescent")
color<-c("blue", "red")
legend(  "topright", legend=Lab,pch=c(22,19),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

moment<-metaMDS(hem.ali, k=3)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.0713509 
## Run 1 stress 0.07135184 
## ... Procrustes: rmse 0.0003176795  max resid 0.0006529927 
## ... Similar to previous best
## Run 2 stress 0.07135107 
## ... Procrustes: rmse 0.0003125567  max resid 0.0006357934 
## ... Similar to previous best
## Run 3 stress 0.09581055 
## Run 4 stress 0.09581243 
## Run 5 stress 0.07135091 
## ... Procrustes: rmse 0.0002325267  max resid 0.0004489137 
## ... Similar to previous best
## Run 6 stress 0.09581176 
## Run 7 stress 0.09581325 
## Run 8 stress 0.07135229 
## ... Procrustes: rmse 0.0005044818  max resid 0.001098918 
## ... Similar to previous best
## Run 9 stress 0.07135072 
## ... New best solution
## ... Procrustes: rmse 0.0002041099  max resid 0.0004146758 
## ... Similar to previous best
## Run 10 stress 0.07135078 
## ... Procrustes: rmse 0.0001761285  max resid 0.0003919795 
## ... Similar to previous best
## Run 11 stress 0.07135076 
## ... Procrustes: rmse 6.081179e-05  max resid 0.000105218 
## ... Similar to previous best
## Run 12 stress 0.07135164 
## ... Procrustes: rmse 0.0003903769  max resid 0.0009003304 
## ... Similar to previous best
## Run 13 stress 0.07135088 
## ... Procrustes: rmse 0.0001376946  max resid 0.0002526987 
## ... Similar to previous best
## Run 14 stress 0.07135159 
## ... Procrustes: rmse 0.000382806  max resid 0.0008377661 
## ... Similar to previous best
## Run 15 stress 0.07135072 
## ... Procrustes: rmse 7.440481e-05  max resid 0.0001546999 
## ... Similar to previous best
## Run 16 stress 0.07135074 
## ... Procrustes: rmse 0.00012075  max resid 0.0002683365 
## ... Similar to previous best
## Run 17 stress 0.0713508 
## ... Procrustes: rmse 5.576722e-05  max resid 9.502438e-05 
## ... Similar to previous best
## Run 18 stress 0.07135075 
## ... Procrustes: rmse 4.769923e-05  max resid 8.960165e-05 
## ... Similar to previous best
## Run 19 stress 0.07135082 
## ... Procrustes: rmse 0.0001276776  max resid 0.0002492422 
## ... Similar to previous best
## Run 20 stress 0.07135201 
## ... Procrustes: rmse 0.0004838851  max resid 0.001138584 
## ... Similar to previous best
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main="Alili Springs")
text(mds.fig, "species",cex=0.5)
points(mds.fig, "sites", pch = 19, col = "green", select = env[env$site=="ALI",]$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="ALI",]$morphotype == "H")


Lab<-c("G", "H")
color<-c("green", "blue")
legend(  "topright", legend=Lab,pch=c(19,15),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

moment<-metaMDS(hem.olaa, k=3)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.06142282 
## Run 1 stress 0.06142317 
## ... Procrustes: rmse 0.0001459225  max resid 0.0002542925 
## ... Similar to previous best
## Run 2 stress 0.06142369 
## ... Procrustes: rmse 0.001569287  max resid 0.002309513 
## ... Similar to previous best
## Run 3 stress 0.06142404 
## ... Procrustes: rmse 0.001889031  max resid 0.003056083 
## ... Similar to previous best
## Run 4 stress 0.06142421 
## ... Procrustes: rmse 0.0004712779  max resid 0.0007472414 
## ... Similar to previous best
## Run 5 stress 0.06142531 
## ... Procrustes: rmse 0.0006532095  max resid 0.0009522755 
## ... Similar to previous best
## Run 6 stress 0.06142253 
## ... New best solution
## ... Procrustes: rmse 0.0004460018  max resid 0.0006143714 
## ... Similar to previous best
## Run 7 stress 0.06177549 
## ... Procrustes: rmse 0.01106618  max resid 0.01887957 
## Run 8 stress 0.06142506 
## ... Procrustes: rmse 0.001144175  max resid 0.001729107 
## ... Similar to previous best
## Run 9 stress 0.06142314 
## ... Procrustes: rmse 0.0004750837  max resid 0.0007324986 
## ... Similar to previous best
## Run 10 stress 0.06142336 
## ... Procrustes: rmse 0.0006656519  max resid 0.0009760464 
## ... Similar to previous best
## Run 11 stress 0.06142309 
## ... Procrustes: rmse 0.0006042221  max resid 0.0008453728 
## ... Similar to previous best
## Run 12 stress 0.06652333 
## Run 13 stress 0.06142247 
## ... New best solution
## ... Procrustes: rmse 0.0003501391  max resid 0.0006710117 
## ... Similar to previous best
## Run 14 stress 0.06652333 
## Run 15 stress 0.06142431 
## ... Procrustes: rmse 0.0008082026  max resid 0.001309491 
## ... Similar to previous best
## Run 16 stress 0.06142235 
## ... New best solution
## ... Procrustes: rmse 0.0001497941  max resid 0.0002665819 
## ... Similar to previous best
## Run 17 stress 0.06652356 
## Run 18 stress 0.06652488 
## Run 19 stress 0.06142374 
## ... Procrustes: rmse 0.0005622079  max resid 0.000919252 
## ... Similar to previous best
## Run 20 stress 0.06142262 
## ... Procrustes: rmse 0.0001982676  max resid 0.0002924331 
## ... Similar to previous best
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main="HAVO Olaa")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 19, col = "green", select = env[env$site=="OLAA",]$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="OLAA",]$morphotype == "H")

Lab<-c("Glabrous", "Hybrid")
color<-c("green", "blue" )
legend(  "topright", legend=Lab,pch=c(19,15),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

# experiment with adding data kaiholena old- only 6 data points
moment<-metaMDS(hem.kaio, k=3)
## Wisconsin double standardization
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.1099498  max resid 0.1686382 
## Run 2 stress 9.310127e-05 
## ... Procrustes: rmse 0.08995409  max resid 0.1557758 
## Run 3 stress 8.792988e-05 
## ... Procrustes: rmse 0.1389952  max resid 0.1699884 
## Run 4 stress 0 
## ... Procrustes: rmse 0.0510628  max resid 0.07329546 
## Run 5 stress 0 
## ... Procrustes: rmse 0.1677124  max resid 0.2726518 
## Run 6 stress 0 
## ... Procrustes: rmse 0.1996001  max resid 0.2891551 
## Run 7 stress 6.345724e-05 
## ... Procrustes: rmse 0.2144208  max resid 0.2853603 
## Run 8 stress 0 
## ... Procrustes: rmse 0.1989755  max resid 0.3117667 
## Run 9 stress 9.697878e-05 
## ... Procrustes: rmse 0.1731629  max resid 0.2866622 
## Run 10 stress 0 
## ... Procrustes: rmse 0.1202304  max resid 0.1895108 
## Run 11 stress 8.569755e-05 
## ... Procrustes: rmse 0.1706865  max resid 0.2845998 
## Run 12 stress 9.006129e-05 
## ... Procrustes: rmse 0.1727758  max resid 0.2864063 
## Run 13 stress 8.904747e-05 
## ... Procrustes: rmse 0.2276194  max resid 0.3055194 
## Run 14 stress 0 
## ... Procrustes: rmse 0.1480634  max resid 0.2556493 
## Run 15 stress 9.791498e-05 
## ... Procrustes: rmse 0.2186505  max resid 0.3068564 
## Run 16 stress 8.368217e-05 
## ... Procrustes: rmse 0.1627448  max resid 0.2784951 
## Run 17 stress 0 
## ... Procrustes: rmse 0.1701545  max resid 0.2513082 
## Run 18 stress 9.282268e-05 
## ... Procrustes: rmse 0.21155  max resid 0.3065631 
## Run 19 stress 0 
## ... Procrustes: rmse 0.1465121  max resid 0.2101937 
## Run 20 stress 0 
## ... Procrustes: rmse 0.1146521  max resid 0.1938358 
## *** No convergence -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(hem.kaio, k = 3): stress is (nearly) zero: you may have
## insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold
mds.fig <- ordiplot(moment, type = "none", main="KAIO limited data")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="KAIO",]$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env[env$site=="KAIO",]$morphotype == "P")
points(mds.fig, "sites", pch = 19, col = "green", select = env[env$site=="KAIO",]$morphotype == "G")

Lab<-c( "Hybrid", "Pubescent", "Glabrous")
color<-c("blue", "red", "green")
legend(  "topright", legend=Lab,pch=c(15,19),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

# explore adding kaio to kaiy. 
hem.kai.all<-rbind(hem.kaiy,hem.kaio)
env.kai.all<-rbind(env[env$site=="KAIY",],env[env$site=="KAIO",])

moment<-metaMDS(hem.kai.all, k=3)
## Wisconsin double standardization
## Run 0 stress 0.1432029 
## Run 1 stress 0.1400188 
## ... New best solution
## ... Procrustes: rmse 0.1367607  max resid 0.328315 
## Run 2 stress 0.1450392 
## Run 3 stress 0.1400853 
## ... Procrustes: rmse 0.1416568  max resid 0.4080269 
## Run 4 stress 0.1349444 
## ... New best solution
## ... Procrustes: rmse 0.1182893  max resid 0.2461261 
## Run 5 stress 0.1349426 
## ... New best solution
## ... Procrustes: rmse 0.0007411143  max resid 0.002142872 
## ... Similar to previous best
## Run 6 stress 0.1369001 
## Run 7 stress 0.1437608 
## Run 8 stress 0.1465123 
## Run 9 stress 0.1401 
## Run 10 stress 0.1450611 
## Run 11 stress 0.1386605 
## Run 12 stress 0.1369041 
## Run 13 stress 0.1349427 
## ... Procrustes: rmse 0.00028015  max resid 0.0005560841 
## ... Similar to previous best
## Run 14 stress 0.1400167 
## Run 15 stress 0.1400186 
## Run 16 stress 0.1465886 
## Run 17 stress 0.1349427 
## ... Procrustes: rmse 0.0006100189  max resid 0.001454125 
## ... Similar to previous best
## Run 18 stress 0.1400879 
## Run 19 stress 0.1400927 
## Run 20 stress 0.1400951 
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main="KAI all")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env.kai.all$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env.kai.all$morphotype == "P")
points(mds.fig, "sites", pch = 19, col = "green", select = env.kai.all$morphotype == "G")

adonis(hem.kai.all~env.kai.all$morphotype, by="margin")
## 
## Call:
## adonis(formula = hem.kai.all ~ env.kai.all$morphotype, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                        Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## env.kai.all$morphotype  2    1.4422 0.72111  2.5311 0.24034  0.001 ***
## Residuals              16    4.5584 0.28490         0.75966           
## Total                  18    6.0006                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(hem.ali~env[env$site=="ALI",]$morphotype, by="margin")
## 
## Call:
## adonis(formula = hem.ali ~ env[env$site == "ALI", ]$morphotype,      by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                     Df SumsOfSqs MeanSqs F.Model      R2
## env[env$site == "ALI", ]$morphotype  1   0.35248 0.35248  1.4706 0.12821
## Residuals                           10   2.39678 0.23968         0.87179
## Total                               11   2.74926                 1.00000
##                                     Pr(>F)
## env[env$site == "ALI", ]$morphotype  0.174
## Residuals                                 
## Total
adonis(hem.kaiy~env[env$site=="KAIY",]$morphotype, by="margin")
## 
## Call:
## adonis(formula = hem.kaiy ~ env[env$site == "KAIY", ]$morphotype,      by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                      Df SumsOfSqs MeanSqs F.Model      R2
## env[env$site == "KAIY", ]$morphotype  1    1.0913 1.09131  4.5097 0.29077
## Residuals                            11    2.6619 0.24199         0.70923
## Total                                12    3.7532                 1.00000
##                                      Pr(>F)   
## env[env$site == "KAIY", ]$morphotype  0.002 **
## Residuals                                     
## Total                                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(hem.ero~env[env$site=="ERO",]$morphotype, by="margin")
## 
## Call:
## adonis(formula = hem.ero ~ env[env$site == "ERO", ]$morphotype,      by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                     Df SumsOfSqs MeanSqs F.Model      R2
## env[env$site == "ERO", ]$morphotype  1    0.4346 0.43459  1.3671 0.12027
## Residuals                           10    3.1789 0.31789         0.87973
## Total                               11    3.6135                 1.00000
##                                     Pr(>F)
## env[env$site == "ERO", ]$morphotype  0.158
## Residuals                                 
## Total
adonis(hem.olaa~env[env$site=="OLAA",]$morphotype, by="margin")
## 
## Call:
## adonis(formula = hem.olaa ~ env[env$site == "OLAA", ]$morphotype,      by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                      Df SumsOfSqs MeanSqs F.Model      R2
## env[env$site == "OLAA", ]$morphotype  1    0.7505 0.75048  2.5677 0.20431
## Residuals                            10    2.9228 0.29228         0.79569
## Total                                11    3.6733                 1.00000
##                                      Pr(>F)  
## env[env$site == "OLAA", ]$morphotype  0.045 *
## Residuals                                    
## Total                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


There are differences in species richness as well as species composition across sites. Diversity is highest at the high productivity sites midway the Big Island chronosequence. Species composition varies between sites, though many species are shared in common.

(ii) Which plant traits account for herbivore community responses to host-plant phenotype?


# partition variance between environment and spatial
           #test<-rda(decostand(trial, "hel"), leaftraits[,c(1:4,6)], xy.2014[,2:3])
#test<-varpart(vegdist(hem), env$z.age, xy[,2:3])
#capscale(vegdist(trial)~ leaftraits[,1:4]+xy[,2:3])
#test
#plot(test)


test<-varpart(vegdist(decostand(hem.2014, "hel")), leaftraits[,c(1,2,3,4,6)], xy.2014[,2:3])
test<-rda(decostand(hem.2014, "hel"), leaftraits[,c(1,2,3,4,6)], xy.2014[,2:3])

test.wout<-rda(decostand(hem.2014, "hel"), leaftraits[,c(1,2,3,4,6)])
test.spat<-rda(decostand(hem.2014, "hel"),  xy.2014[,2:3])
anova(test)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(X = decostand(hem.2014, "hel"), Y = leaftraits[, c(1, 2, 3, 4, 6)], Z = xy.2014[, 2:3])
##          Df Variance      F Pr(>F)   
## Model     6  0.22081 1.7192  0.008 **
## Residual 15  0.32109                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(hem.2014~leaftraits$percN+leaftraits$meanWatercontent+leaftraits$meanSLA+leaftraits$percP+xy.2014$long+env.2014$morphotype+xy.2014$lat, by="margin")
## 
## Call:
## adonis(formula = hem.2014 ~ leaftraits$percN + leaftraits$meanWatercontent +      leaftraits$meanSLA + leaftraits$percP + xy.2014$long + env.2014$morphotype +      xy.2014$lat, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                             Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## leaftraits$percN             1    1.1225 1.12250  4.4231 0.13512  0.001
## leaftraits$meanWatercontent  1    0.4434 0.44340  1.7472 0.05338  0.054
## leaftraits$meanSLA           1    0.3591 0.35912  1.4151 0.04323  0.161
## leaftraits$percP             1    0.9557 0.95568  3.7658 0.11504  0.001
## xy.2014$long                 1    0.5235 0.52347  2.0627 0.06301  0.025
## env.2014$morphotype          2    0.6836 0.34179  1.3468 0.08229  0.136
## xy.2014$lat                  1    0.4127 0.41273  1.6263 0.04968  0.086
## Residuals                   15    3.8067 0.25378         0.45824       
## Total                       23    8.3072                 1.00000       
##                                
## leaftraits$percN            ***
## leaftraits$meanWatercontent .  
## leaftraits$meanSLA             
## leaftraits$percP            ***
## xy.2014$long                *  
## env.2014$morphotype            
## xy.2014$lat                 .  
## Residuals                      
## Total                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(hem.2014~xy.2014$lat+xy.2014$long+leaftraits$percN+leaftraits$percP+leaftraits$meanSLA+leaftraits$meanWatercontent+env.2014$morphotype, by="margin")
## 
## Call:
## adonis(formula = hem.2014 ~ xy.2014$lat + xy.2014$long + leaftraits$percN +      leaftraits$percP + leaftraits$meanSLA + leaftraits$meanWatercontent +      env.2014$morphotype, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                             Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## xy.2014$lat                  1    1.1255 1.12547  4.4349 0.13548  0.001
## xy.2014$long                 1    0.7394 0.73945  2.9137 0.08901  0.005
## leaftraits$percN             1    0.9300 0.93003  3.6647 0.11196  0.001
## leaftraits$percP             1    0.3141 0.31412  1.2378 0.03781  0.275
## leaftraits$meanSLA           1    0.3360 0.33602  1.3241 0.04045  0.195
## leaftraits$meanWatercontent  1    0.4945 0.49448  1.9485 0.05952  0.034
## env.2014$morphotype          2    0.5609 0.28045  1.1051 0.06752  0.341
## Residuals                   15    3.8067 0.25378         0.45824       
## Total                       23    8.3072                 1.00000       
##                                
## xy.2014$lat                 ***
## xy.2014$long                ** 
## leaftraits$percN            ***
## leaftraits$percP               
## leaftraits$meanSLA             
## leaftraits$meanWatercontent *  
## env.2014$morphotype            
## Residuals                      
## Total                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# rda/varpart: when testing only 2014 hybrid zone sites from 2014, xy coordinates explain some. better explained still by leaftraits
# adonis seems to indicate latitude is important predictor- sig even after all other leaftrait variables
# while almost always sig (for 170 data points and 24), R2 is less than 0.1 in all cases  


adonis(hem.2014~leaftraits$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$meanSLA+leaftraits$percP+scale(xy.2014$lat,scale=T,center=T), by="margin")
## 
## Call:
## adonis(formula = hem.2014 ~ leaftraits$morphotype + leaftraits$percN +      leaftraits$meanWatercontent + leaftraits$meanSLA + leaftraits$percP +      scale(xy.2014$lat, scale = T, center = T), by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                           Df SumsOfSqs MeanSqs F.Model
## leaftraits$morphotype                      2    1.3702 0.68512  2.6179
## leaftraits$percN                           1    0.6270 0.62697  2.3957
## leaftraits$meanWatercontent                1    0.5596 0.55960  2.1383
## leaftraits$meanSLA                         1    0.3781 0.37811  1.4448
## leaftraits$percP                           1    0.6230 0.62301  2.3806
## scale(xy.2014$lat, scale = T, center = T)  1    0.5619 0.56193  2.1472
## Residuals                                 16    4.1873 0.26171        
## Total                                     23    8.3072                
##                                                R2 Pr(>F)    
## leaftraits$morphotype                     0.16495  0.001 ***
## leaftraits$percN                          0.07547  0.006 ** 
## leaftraits$meanWatercontent               0.06736  0.032 *  
## leaftraits$meanSLA                        0.04552  0.153    
## leaftraits$percP                          0.07500  0.012 *  
## scale(xy.2014$lat, scale = T, center = T) 0.06764  0.036 *  
## Residuals                                 0.50406           
## Total                                     1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(hem.2014~env.2014$tree.height+leaftraits$morphotype+leaftraits$meanSLA+leaftraits$percN+env.2014$z.age+xy.2014$lat+xy.2014$long+env.2014$site, by="margin")
## 
## Call:
## adonis(formula = hem.2014 ~ env.2014$tree.height + leaftraits$morphotype +      leaftraits$meanSLA + leaftraits$percN + env.2014$z.age +      xy.2014$lat + xy.2014$long + env.2014$site, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                       Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## env.2014$tree.height   1    0.6339 0.63394 2.37893 0.07631  0.012 * 
## leaftraits$morphotype  2    1.0796 0.53981 2.02572 0.12996  0.014 * 
## leaftraits$meanSLA     1    0.7449 0.74495 2.79552 0.08968  0.005 **
## leaftraits$percN       1    0.3277 0.32774 1.22988 0.03945  0.251   
## env.2014$z.age         1    0.4335 0.43353 1.62688 0.05219  0.099 . 
## xy.2014$lat            1    0.7524 0.75241 2.82353 0.09057  0.006 **
## xy.2014$long           1    0.3657 0.36574 1.37249 0.04403  0.208   
## env.2014$site          2    0.5050 0.25250 0.94754 0.06079  0.542   
## Residuals             13    3.4642 0.26648         0.41702          
## Total                 23    8.3072                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### try again for all 2014, not just hybrid

env.sub<-env[env$vialID%in%rownames(leaftraits.full),]
hem.sub<-hem[rownames(hem)%in%rownames(leaftraits.full),]

adonis(hem.sub~leaftraits.full$morphotype+leaftraits.full$meanSLA+leaftraits.full$percN+env.sub$z.age+env.sub$site, by="margin")
## 
## Call:
## adonis(formula = hem.sub ~ leaftraits.full$morphotype + leaftraits.full$meanSLA +      leaftraits.full$percN + env.sub$z.age + env.sub$site, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## leaftraits.full$morphotype  2    2.0183 1.00917  3.7114 0.07418  0.001 ***
## leaftraits.full$meanSLA     1    1.2077 1.20774  4.4417 0.04439  0.001 ***
## leaftraits.full$percN       1    0.2891 0.28909  1.0632 0.01062  0.369    
## env.sub$z.age               1    0.8705 0.87046  3.2013 0.03199  0.001 ***
## env.sub$site               11    7.5954 0.69049  2.5394 0.27916  0.001 ***
## Residuals                  56   15.2271 0.27191         0.55965           
## Total                      72   27.2081                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
temp<-hem[env$sampleyear=="2014",]
temp<-temp[env$vialID%in%rownames(leaftraits.full),]

mod.leaf<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+leaftraits$meanWatercontent+leaftraits$percP)
mod.leaf.age<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$age)
mod.leaf.tree<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height)
mod.leaf.age.site<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
mod.leaf.site<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)


modfull<-cca(hem[rownames(hem)%in%rownames(leaftraits),]~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m1<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m2<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+leaftraits$percP+leaftraits$meanWatercontent)
m3<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m4<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)


m5<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$site)
m6<-cca(hem.2014~leaftraits$morphotype+leaftraits$meanSLA+env.2014$site )       
  anova(m5, by = "margin", step=200 )      
## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = hem.2014 ~ leaftraits$percN + leaftraits$morphotype + leaftraits$meanSLA + env.2014$site)
##                       Df ChiSquare      F Pr(>F)   
## leaftraits$percN       1   0.09340 1.0262  0.365   
## leaftraits$morphotype  2   0.25951 1.4257  0.044 * 
## leaftraits$meanSLA     1   0.11201 1.2307  0.204   
## env.2014$site          3   0.48967 1.7934  0.004 **
## Residual              16   1.45618                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  anova(m5,m6)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: hem.2014 ~ leaftraits$percN + leaftraits$morphotype + leaftraits$meanSLA + env.2014$site
## Model 2: hem.2014 ~ leaftraits$morphotype + leaftraits$meanSLA + env.2014$site
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)
## 1    16       1.4562                           
## 2    17       1.5496 -1 -0.093399 1.0262   0.36
m7<-cca(hem.2014~leaftraits$morphotype+env.2014$site )   
anova(m6,m7)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: hem.2014 ~ leaftraits$morphotype + leaftraits$meanSLA + env.2014$site
## Model 2: hem.2014 ~ leaftraits$morphotype + env.2014$site
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)  
## 1    17       1.5496                             
## 2    18       1.6780 -1  -0.12845 1.4091  0.092 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m8<-cca(hem.2014~env.2014$site )  
anova(m7,m8)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: hem.2014 ~ leaftraits$morphotype + env.2014$site
## Model 2: hem.2014 ~ env.2014$site
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)    
## 1    18       1.6780                               
## 2    20       2.1829 -2  -0.50488 2.7079  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model with just site is best

Do CCA and permanova with data lumped for 2014 and 2015

# subset data in species by site matrix such that it has only 24 rows, one for each site plot combination
# each row should be the sum of collecting in 2014 and 2015

# be careful to exclude anywhere where sampling effort is wildly different

hem$vialID<-rownames(hem)
temp<-melt(hem)
## Using vialID as id variables
#trial<-dcast(temp, site+plot~variable, fun.aggregate = sum,value.var="value")

trial.env<-merge(temp, env[,1:3], by="vialID")
trial<-dcast(trial.env, site+plot~variable, fun.aggregate = sum,value.var="value")

trial<-trial[trial$site%in%env.both$site,]

hem<-hem[,-84]
rowSums(hem[env$site=="KAIY",])
##  25  26  28  30  42  58 141 142 143 163 164 165 166 
##   6  38  12   4   4  10   8  22  16  76  18  29  25
rowSums(trial[,-c(1:2)])
##   1   2   3   4   5   6   7   8   9  10  11  12  25  26  27  28  29  30 
##  42  36 130  83  31  12  44  42  65  38  20  18  22  22  12  39  37 136 
##  73  74  75  76  77  78 
##  55 611  56 592 120 494
##################################################

env.trial<-trial[,1:2]
env.trial<-unique(merge(env.trial, env[,c(2,3,10)], by=c("site","plot")))
trial<-trial[,-c(1:2)]


moment<-metaMDS(trial, k=3)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1492291 
## Run 1 stress 0.1659777 
## Run 2 stress 0.1489775 
## ... New best solution
## ... Procrustes: rmse 0.03255228  max resid 0.1093367 
## Run 3 stress 0.1661723 
## Run 4 stress 0.1489705 
## ... New best solution
## ... Procrustes: rmse 0.06705354  max resid 0.1853576 
## Run 5 stress 0.148974 
## ... Procrustes: rmse 0.0009068024  max resid 0.002178237 
## ... Similar to previous best
## Run 6 stress 0.1489704 
## ... New best solution
## ... Procrustes: rmse 0.0012843  max resid 0.002781557 
## ... Similar to previous best
## Run 7 stress 0.1489776 
## ... Procrustes: rmse 0.06670506  max resid 0.1817755 
## Run 8 stress 0.1489698 
## ... New best solution
## ... Procrustes: rmse 0.001356732  max resid 0.004352775 
## ... Similar to previous best
## Run 9 stress 0.1489692 
## ... New best solution
## ... Procrustes: rmse 0.001615042  max resid 0.004077342 
## ... Similar to previous best
## Run 10 stress 0.1690282 
## Run 11 stress 0.1800283 
## Run 12 stress 0.1492302 
## ... Procrustes: rmse 0.05076003  max resid 0.1856264 
## Run 13 stress 0.1489787 
## ... Procrustes: rmse 0.06748909  max resid 0.1828223 
## Run 14 stress 0.16617 
## Run 15 stress 0.1489843 
## ... Procrustes: rmse 0.0672218  max resid 0.1824667 
## Run 16 stress 0.1489705 
## ... Procrustes: rmse 0.0004047275  max resid 0.00093899 
## ... Similar to previous best
## Run 17 stress 0.1661743 
## Run 18 stress 0.1489707 
## ... Procrustes: rmse 0.00210368  max resid 0.004964002 
## ... Similar to previous best
## Run 19 stress 0.1492229 
## ... Procrustes: rmse 0.05149002  max resid 0.185743 
## Run 20 stress 0.1489796 
## ... Procrustes: rmse 0.06739322  max resid 0.1828275 
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main=" all")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env.trial$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env.trial$morphotype == "P")
points(mds.fig, "sites", pch = 19, col = "green", select = env.trial$morphotype == "G")

adonis(trial~env.trial$morphotype, by="margin")
## 
## Call:
## adonis(formula = trial ~ env.trial$morphotype, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## env.trial$morphotype  2    1.3344 0.66719  2.4863 0.19146  0.001 ***
## Residuals            21    5.6352 0.26834         0.80854           
## Total                23    6.9696                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# with 2014 and 2015 data summed pretty similar results to when with twice as many data points

# pretend it has vialID
env.trial<-merge(env.trial, env.2014[,1:3], by=c("site","plot"))

rownames(trial)<-env.trial$vialID

# make sure all same order
env.trial<-env.trial[order(env.trial$vialID),]
trial<-trial[order(rownames(trial)),]


m1<-adonis(trial~leaftraits$meanWatercontent+env.trial$morphotype+leaftraits$meanSLA+leaftraits$percN+leaftraits$percP, by="margin")
# even as first variable water content is not sig
adonis(trial~env.trial$morphotype+leaftraits$meanSLA+leaftraits$percN+leaftraits$percP, by="margin")
## 
## Call:
## adonis(formula = trial ~ env.trial$morphotype + leaftraits$meanSLA +      leaftraits$percN + leaftraits$percP, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## env.trial$morphotype  2    1.3369 0.66845 2.73289 0.19182  0.001 ***
## leaftraits$meanSLA    1    0.5869 0.58685 2.39928 0.08420  0.014 *  
## leaftraits$percN      1    0.2387 0.23869 0.97584 0.03425  0.470    
## leaftraits$percP      1    0.4044 0.40440 1.65336 0.05802  0.085 .  
## Residuals            18    4.4027 0.24460         0.63171           
## Total                23    6.9696                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# percN NS in third place, what about up front
adonis(trial~leaftraits$percN+leaftraits$meanSLA+env.trial$morphotype+leaftraits$percP, by="margin")
## 
## Call:
## adonis(formula = trial ~ leaftraits$percN + leaftraits$meanSLA +      env.trial$morphotype + leaftraits$percP, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## leaftraits$percN      1    0.7208 0.72078  2.9468 0.10342  0.002 **
## leaftraits$meanSLA    1    0.3543 0.35427  1.4484 0.05083  0.156   
## env.trial$morphotype  2    1.0874 0.54370  2.2229 0.15602  0.004 **
## leaftraits$percP      1    0.4044 0.40440  1.6534 0.05802  0.092 . 
## Residuals            18    4.4027 0.24460         0.63171          
## Total                23    6.9696                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# does better in first place, seems to correlate alot with SLA. what if only including SLA- sig in last place?
adonis(trial~env.trial$morphotype+leaftraits$percP+leaftraits$meanSLA, by="margin")
## 
## Call:
## adonis(formula = trial ~ env.trial$morphotype + leaftraits$percP +      leaftraits$meanSLA, by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## env.trial$morphotype  2    1.3369 0.66845  2.6954 0.19182  0.002 **
## leaftraits$percP      1    0.2545 0.25451  1.0263 0.03652  0.425   
## leaftraits$meanSLA    1    0.6662 0.66618  2.6862 0.09558  0.004 **
## Residuals            19    4.7120 0.24800         0.67608          
## Total                23    6.9696                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# SLA still sig. percP is not in second place (it is sig in first place)
adonis(trial~leaftraits$meanSLA+env.trial$morphotype, by="margin")
## 
## Call:
## adonis(formula = trial ~ leaftraits$meanSLA + env.trial$morphotype,      by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## leaftraits$meanSLA    1    0.8404 0.84042  3.3311 0.12058  0.002 **
## env.trial$morphotype  2    1.0833 0.54167  2.1470 0.15544  0.013 * 
## Residuals            20    5.0458 0.25229         0.72398          
## Total                23    6.9696                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mean SLA and morphotype sig
# regardless of which comes first, together R2 0.28, in both cases ~p<0.01

# what about xy:
adonis(trial~leaftraits$percN+leaftraits$percP+env.trial$morphotype+leaftraits$meanSLA+xy.2014$lat, by="margin")
## 
## Call:
## adonis(formula = trial ~ leaftraits$percN + leaftraits$percP +      env.trial$morphotype + leaftraits$meanSLA + xy.2014$lat,      by = "margin") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## leaftraits$percN      1    0.7208 0.72078  3.1357 0.10342  0.002 **
## leaftraits$percP      1    0.5768 0.57683  2.5095 0.08276  0.010 **
## env.trial$morphotype  2    0.8222 0.41112  1.7886 0.11798  0.037 * 
## leaftraits$meanSLA    1    0.4470 0.44700  1.9446 0.06414  0.024 * 
## xy.2014$lat           1    0.4951 0.49509  2.1539 0.07104  0.019 * 
## Residuals            17    3.9076 0.22986         0.56067          
## Total                23    6.9696                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod4<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$lat+xy.2014$long)
mod3<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent)
mod5<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$lat)
anova(mod3,mod4,mod5)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 3: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)  
## 1    18       1.9822                             
## 2    15       1.5024  3   0.47982 1.5968  0.015 *
## 3    16       1.6095 -1  -0.10710 1.0693  0.377  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod3,mod4)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)  
## 1    18       1.9822                             
## 2    15       1.5024  3   0.47982 1.5968  0.012 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model shortest is best
# try removing latitude
mod6<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$long)
anova(mod4,mod6)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$long
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)
## 1    15       1.5024                           
## 2    16       1.5896 -1 -0.087169 0.8703  0.642
# shortest model is best
mod3<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent)
mod2<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN)
anova(mod2,mod3)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)   
## 1    19       2.2269                              
## 2    18       1.9822  1   0.24468 2.2219  0.007 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model 3 seems better
mod<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$meanWatercontent)
anova(mod,mod3)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)
## 1    19       2.1035                           
## 2    18       1.9822  1   0.12131 1.1016  0.357
# no difference
mod6<-cca(decostand(trial, "total")~leaftraits$meanSLA+leaftraits$meanWatercontent)
anova(mod,mod6)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + leaftraits$meanWatercontent
##   ResDf ResChiSquare Df ChiSquare     F Pr(>F)   
## 1    19       2.1035                             
## 2    21       2.5572 -2  -0.45369 2.049  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod7<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype)
anova(mod,mod7)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)   
## 1    19       2.1035                              
## 2    20       2.3389 -1  -0.23535 2.1258  0.006 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod8<-cca(decostand(trial, "total")~leaftraits$meanSLA)
anova(mod4,mod7,mod8)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype
## Model 3: decostand(trial, "total") ~ leaftraits$meanSLA
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)    
## 1    15       1.5024                               
## 2    20       2.3389 -5  -0.83648 1.6703  0.002 ** 
## 3    22       2.7205 -2  -0.38163 1.9051  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod9<-cca(decostand(trial, "total")~env.trial$morphotype)
anova(mod4,mod7,mod9)
## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype
## Model 3: decostand(trial, "total") ~ env.trial$morphotype
##   ResDf ResChiSquare Df ChiSquare      F Pr(>F)    
## 1    15       1.5024                               
## 2    20       2.3389 -5  -0.83648 1.6703  0.001 ***
## 3    21       2.5950 -1  -0.25611 2.5571  0.003 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# can't compare model with just morphotype or just SLA, but F value in the same comparison (full model, model with the two variables and model with just one) is higher for morphotype
# best model mod9




mod4<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$lat+xy.2014$long)

plot(mod4, display = "sites", scaling = 3, type="points")
text(mod4, scaling = 3, display = "bp", pch=4) 
text(mod4, display = "species", scaling = 3, col="red", cex=0.5)

anova(mod4, by = "margin", step=200 )
## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long)
##                             Df ChiSquare      F Pr(>F)  
## leaftraits$meanSLA           1   0.14319 1.4297  0.146  
## env.trial$morphotype         2   0.30895 1.5423  0.037 *
## leaftraits$percN             1   0.11347 1.1329  0.299  
## leaftraits$meanWatercontent  1   0.14256 1.4233  0.124  
## leaftraits$percP             1   0.11941 1.1922  0.303  
## xy.2014$lat                  1   0.08717 0.8703  0.633  
## xy.2014$long                 1   0.10710 1.0693  0.397  
## Residual                    15   1.50239                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod4
## Call: cca(formula = decostand(trial, "total") ~ leaftraits$meanSLA
## + env.trial$morphotype + leaftraits$percN +
## leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat +
## xy.2014$long)
## 
##               Inertia Proportion Rank
## Total          3.0213     1.0000     
## Constrained    1.5189     0.5027    8
## Unconstrained  1.5024     0.4973   15
## Inertia is scaled Chi-square 
## 38 species (variables) deleted due to missingness
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8 
## 0.4530 0.3274 0.2532 0.2008 0.1061 0.0760 0.0602 0.0422 
## 
## Eigenvalues for unconstrained axes:
##     CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8     CA9 
## 0.31373 0.27569 0.16823 0.16084 0.13136 0.09011 0.08342 0.06841 0.05593 
##    CA10    CA11    CA12    CA13    CA14    CA15 
## 0.04687 0.03611 0.02644 0.01971 0.01593 0.00962
# morphotype, mean water and percP are sig, SLA marginally
# on figure, percN, percP and SLA completely overlaying eachother. morphotype H and P also fall in same spot
# explains 41% variation



# explore:
# seems like all of the dimensions around leaftraits are captured by meanSLA (cor strongly with percN, percP), morphotype and water content
vare.cca <- cca(decostand(trial, "total") ~ env.trial$morphotype) 

plot(vare.cca, display = "sites", scaling = 3, type="points")
text(vare.cca, scaling = 3, display = "bp", pch=4) 
text(vare.cca, display = "species", scaling = 3, col="red", cex=0.5)

anova(vare.cca, by = "margin", step=200 )
## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = decostand(trial, "total") ~ env.trial$morphotype)
##                      Df ChiSquare     F Pr(>F)    
## env.trial$morphotype  2   0.42632 1.725  0.001 ***
## Residual             21   2.59498                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vare.cca
## Call: cca(formula = decostand(trial, "total") ~
## env.trial$morphotype)
## 
##               Inertia Proportion Rank
## Total          3.0213     1.0000     
## Constrained    0.4263     0.1411    2
## Unconstrained  2.5950     0.8589   21
## Inertia is scaled Chi-square 
## 38 species (variables) deleted due to missingness
## 
## Eigenvalues for constrained axes:
##    CCA1    CCA2 
## 0.23183 0.19449 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.5219 0.3710 0.2926 0.2466 0.2218 0.1693 0.1407 0.1262 
## (Showed only 8 of all 21 unconstrained eigenvalues)
# together 30% variation


(iii) Which insect traits account for herbivore community responses to host-plant phenotype



Are communities distinctly different between different morphotype trees across all youngish sites?

moment<-metaMDS(hem.both, k=3)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1790701 
## Run 1 stress 0.1860155 
## Run 2 stress 0.1839765 
## Run 3 stress 0.1892048 
## Run 4 stress 0.1865272 
## Run 5 stress 0.1791475 
## ... Procrustes: rmse 0.008537058  max resid 0.04929761 
## Run 6 stress 0.1832359 
## Run 7 stress 0.1791218 
## ... Procrustes: rmse 0.01407452  max resid 0.08428514 
## Run 8 stress 0.1889459 
## Run 9 stress 0.179072 
## ... Procrustes: rmse 0.0007812932  max resid 0.002826476 
## ... Similar to previous best
## Run 10 stress 0.1790731 
## ... Procrustes: rmse 0.0008433542  max resid 0.00345565 
## ... Similar to previous best
## Run 11 stress 0.1832108 
## Run 12 stress 0.1826911 
## Run 13 stress 0.1860141 
## Run 14 stress 0.1854629 
## Run 15 stress 0.186774 
## Run 16 stress 0.1790705 
## ... Procrustes: rmse 0.0008900677  max resid 0.004086788 
## ... Similar to previous best
## Run 17 stress 0.1791349 
## ... Procrustes: rmse 0.01539417  max resid 0.08989465 
## Run 18 stress 0.1805679 
## Run 19 stress 0.1791347 
## ... Procrustes: rmse 0.01156497  max resid 0.06346796 
## Run 20 stress 0.1791222 
## ... Procrustes: rmse 0.01373238  max resid 0.08235413 
## *** Solution reached
mds.fig <- ordiplot(moment, type = "none", main="all hybrid zone sites")
#text(mds.fig, "species")
points(mds.fig, "sites", pch = 19, col = "green", select = env.both$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env.both$morphotype == "H")
points(mds.fig, "sites", pch = 17, col = "red", select = env.both$morphotype == "P")

Lab<-c("Glabrous", "Hybrid", "Pubescent")
color<-c("green", "blue", "red")
legend(  "topright", legend=Lab,pch=c(19,15,17),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

## across all sites (n=170)
moment<-metaMDS(hem, k=3)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2143791 
## Run 1 stress 0.2144895 
## ... Procrustes: rmse 0.02158443  max resid 0.09391575 
## Run 2 stress 0.2143115 
## ... New best solution
## ... Procrustes: rmse 0.04348152  max resid 0.1854511 
## Run 3 stress 0.2143784 
## ... Procrustes: rmse 0.03803577  max resid 0.1705378 
## Run 4 stress 0.2147627 
## ... Procrustes: rmse 0.02768923  max resid 0.1083661 
## Run 5 stress 0.2148255 
## Run 6 stress 0.2135422 
## ... New best solution
## ... Procrustes: rmse 0.03912804  max resid 0.1809095 
## Run 7 stress 0.2141577 
## Run 8 stress 0.2132945 
## ... New best solution
## ... Procrustes: rmse 0.01515823  max resid 0.08904171 
## Run 9 stress 0.2163937 
## Run 10 stress 0.2144438 
## Run 11 stress 0.2133151 
## ... Procrustes: rmse 0.002346661  max resid 0.0125573 
## Run 12 stress 0.2142829 
## Run 13 stress 0.2134329 
## ... Procrustes: rmse 0.01233698  max resid 0.1037661 
## Run 14 stress 0.2135966 
## ... Procrustes: rmse 0.01799393  max resid 0.08714684 
## Run 15 stress 0.2148402 
## Run 16 stress 0.2134204 
## ... Procrustes: rmse 0.01341335  max resid 0.07907385 
## Run 17 stress 0.2146835 
## Run 18 stress 0.2144819 
## Run 19 stress 0.2180064 
## Run 20 stress 0.2142619 
## *** No convergence -- monoMDS stopping criteria:
##      6: no. of iterations >= maxit
##     14: stress ratio > sratmax
mds.fig <- ordiplot(moment, type = "none", main="all sites")
text(mds.fig, "species",cex=0.5)
points(mds.fig, "sites", pch = 19, col = "grey", select = env$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env$morphotype == "H")
points(mds.fig, "sites", pch = 17, col = "red", select = env$morphotype == "P")

Lab<-c("Glabrous", "Hybrid", "Pubescent")
color<-c("grey", "blue", "red")
legend(  "topright", legend=Lab,pch=c(19,15,17),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)

# for all sites without double measures (n=97)
#moment<-metaMDS(hem.all, k=3)
#mds.fig <- ordiplot(moment, type = "none", main="all sites")
#text(mds.fig, "species",cex=0.5)
#points(mds.fig, "sites", pch = 19, col = "grey", select = env$morphotype == "G")
#points(mds.fig, "sites", pch = 15, col = "blue", select = env$morphotype == "H")
#points(mds.fig, "sites", pch = 17, col = "red", select = env$morphotype == "P")


phosphorus and morphotype are only sig when only/first term

mean SLA and mean water content are always sig no matter their placement

perc N significant unless after SLA and before water content


What does diversity and abundance look like across these sites with multiple phenotypes?

plot(1/(1-ChaoSimpson(as.data.frame(t(hem.ero)) )[,2])~env[env$site=="ERO",]$morphotype)

plot(1/(1-ChaoSimpson(as.data.frame(t(hem.kaiy)) )[,2])~env[env$site=="KAIY",]$morphotype)
## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (Inf) in boxplot 4 is not drawn

plot(1/(1-ChaoSimpson(as.data.frame(t(hem.olaa)) )[,2])~env[env$site=="OLAA",]$morphotype)

plot(1/(1-ChaoSimpson(as.data.frame(t(hem.ali)) )[,2])~env[env$site=="ALI",]$morphotype)
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

plot(ChaoRichness(as.data.frame(t(hem.ero)) )[,2]~env[env$site=="ERO",]$morphotype)

plot(ChaoRichness(as.data.frame(t(hem.kaiy)) )[,2]~env[env$site=="KAIY",]$morphotype)

plot(ChaoRichness(as.data.frame(t(hem.olaa)) )[,2]~env[env$site=="OLAA",]$morphotype)

plot(ChaoRichness(as.data.frame(t(hem.ali)) )[,2]~env[env$site=="ALI",]$morphotype)

plot(rowSums(hem.ero)~env[env$site=="ERO",]$morphotype)

plot(rowSums(hem.kaiy)~env[env$site=="KAIY",]$morphotype)

plot(rowSums(hem.olaa)~env[env$site=="OLAA",]$morphotype)

plot(rowSums(hem.ali)~env[env$site=="ALI",]$morphotype)

anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.ero)) )[,2])~env[env$site=="ERO",]$morphotype))
## Analysis of Variance Table
## 
## Response: 1/(1 - ChaoSimpson(as.data.frame(t(hem.ero)))[, 2])
##                                     Df  Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "ERO", ]$morphotype  1   9.503  9.5031  0.9171 0.3608
## Residuals                           10 103.626 10.3626
#anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.kaiy)) )[,2])~env[env$site=="KAIY",]$morphotype))
anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.olaa)) )[,2])~env[env$site=="OLAA",]$morphotype))
## Analysis of Variance Table
## 
## Response: 1/(1 - ChaoSimpson(as.data.frame(t(hem.olaa)))[, 2])
##                                      Df  Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "OLAA", ]$morphotype  1  0.2931  0.2931  0.1316 0.7244
## Residuals                            10 22.2796  2.2280
anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.ali)) )[,2])~env[env$site=="ALI",]$morphotype))
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Analysis of Variance Table
## 
## Response: 1/(1 - ChaoSimpson(as.data.frame(t(hem.ali)))[, 2])
##                                     Df  Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "ALI", ]$morphotype  1  46.116  46.116  1.8435 0.2044
## Residuals                           10 250.158  25.016
anova(lm(rowSums(hem.kaiy)~env[env$site=="KAIY",]$morphotype))
## Analysis of Variance Table
## 
## Response: rowSums(hem.kaiy)
##                                      Df Sum Sq Mean Sq F value  Pr(>F)  
## env[env$site == "KAIY", ]$morphotype  1 1589.0 1589.00   5.803 0.03468 *
## Residuals                            11 3012.1  273.83                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p<0.05
anova(lm(rowSums(hem.olaa)~env[env$site=="OLAA",]$morphotype))
## Analysis of Variance Table
## 
## Response: rowSums(hem.olaa)
##                                      Df Sum Sq Mean Sq F value  Pr(>F)  
## env[env$site == "OLAA", ]$morphotype  1  81667   81667  4.6195 0.05713 .
## Residuals                            10 176786   17679                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p<0.1


# Across all 12 sites:
env.both$morphotype<-factor(env.both$morphotype, levels=c("H", "G", "P"))
plot(rowSums(hem.both)~env.both$morphotype, ylim=c(0,300), ylab="Abundance")

anova(lm(rowSums(hem.both)~env.both$morphotype))
## Analysis of Variance Table
## 
## Response: rowSums(hem.both)
##                     Df Sum Sq Mean Sq F value   Pr(>F)   
## env.both$morphotype  2  91547   45773  5.9593 0.004996 **
## Residuals           46 353328    7681                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(rowSums(hem.both)~env.both$morphotype))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = rowSums(hem.both) ~ env.both$morphotype)
## 
## $`env.both$morphotype`
##           diff         lwr       upr     p adj
## G-H  82.796154    7.178384 158.41392 0.0289491
## P-H  -8.841346  -88.095293  70.41260 0.9605958
## P-G -91.637500 -162.829276 -20.44572 0.0086619
# p<0.05

tem<-1-ChaoSimpson(as.data.frame(t(hem.both), na.rm=TRUE) )
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
#tem[10,2]<-1
#tem[65,2]<-1
tem<-ifelse(tem[,2]==0,0,1/tem[,2])
anova(lm(tem~env.both$morphotype))
## Analysis of Variance Table
## 
## Response: tem
##                     Df Sum Sq Mean Sq F value  Pr(>F)  
## env.both$morphotype  2  68.32  34.160  2.6805 0.07922 .
## Residuals           46 586.22  12.744                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmer(tem~env.both$morphotype+(1|env.both$site)))
## Analysis of Variance Table
##                     Df Sum Sq Mean Sq F value
## env.both$morphotype  2 53.867  26.933  2.2435
tem<-exp(ChaoShannon(as.data.frame(t(hem.both)) )[,2])
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
#tem[2]<-1
#tem[26]<-1 
anova(lm(tem~env.both$morphotype))
## Analysis of Variance Table
## 
## Response: tem
##                     Df Sum Sq Mean Sq F value  Pr(>F)  
## env.both$morphotype  2  68.43  34.217   3.074 0.05584 .
## Residuals           46 512.04  11.131                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmer(tem~env.both$morphotype+(1|env.both$site)))
## Analysis of Variance Table
##                     Df Sum Sq Mean Sq F value
## env.both$morphotype  2  66.67  33.335  3.0035
tem<-ChaoRichness(as.data.frame(t(hem.both)))[,2]
anova(lm(tem~env.both$morphotype))
## Analysis of Variance Table
## 
## Response: tem
##                     Df  Sum Sq Mean Sq F value  Pr(>F)  
## env.both$morphotype  2  233.82 116.912  2.6616 0.08058 .
## Residuals           46 2020.60  43.926                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmer(tem~env.both$morphotype+(1|env.both$site)))
## Analysis of Variance Table
##                     Df Sum Sq Mean Sq F value
## env.both$morphotype  2 241.78  120.89  2.7891
plot(tem~env.both$morphotype, ylab="Species Richness")

# richness does diff sig between morphotype, diversity shannon or simpson does not.
# ie the difference between glabrous and hybrid and glabrous and pubescent is in rare species- glabrous has more rare species than the other two morphotypes, but once you weigh rare species less, there are no differences in diversity



### plots for all sites
plot(diversity(hem)/log(specnumber(hem))~env$morphotype, ylab="Pielou's Evenness")

plot(rowSums(hem)~env$morphotype,  ylab="Abundance",xlab="Tree morphotype", main="across big island (n=170)", ylim=c(0,350))

tem<-exp(ChaoShannon(as.data.frame(t(hem)) )[,2])
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.

## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
plot(tem~env$morphotype, ylab="Exp Shannon Index",main="across big island (n=170)", xlab="Tree morphotype")


How much variation is explained by spatial proximity?


GLMER
# dataset with 2014 and 2015 data
hem.both$vialID<-rownames(hem.both)
hem.both.long<-melt(hem.both)
## Using vialID as id variables
hem.both.long<-merge(env.both, hem.both.long, by="vialID")
leaftraits$vialID<-rownames(leaftraits)
temp<-merge(leaftraits,env[,c(1:3)], by="vialID")
temp<-temp[,-1]   
hem.both.long<-hem.both.long[,-1]

hem.both.long<-merge(hem.both.long, temp, by=c("site", "plot"))

hem.both.long$percN.standard<-(hem.both.long$percN-mean(hem.both.long$percN))/sd(hem.both.long$percN)
hem.both.long$percP.standard<-(hem.both.long$percP-mean(hem.both.long$percP))/sd(hem.both.long$percP)
hem.both.long$meanSLA.standard<-(hem.both.long$meanSLA-mean(hem.both.long$meanSLA))/sd(hem.both.long$meanSLA)
hem.both.long$meanWatercontent.standard<-(hem.both.long$meanWatercontent-mean(hem.both.long$meanWatercontent))/sd(hem.both.long$meanWatercontent)


############## analysis data 2014, 2015

modelP=glmer(value ~morphotype.x +meanSLA.standard +percP.standard+(1+percP.standard|variable)+(1|site),  data    =hem.both.long,   family    =poisson)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00111975
## (tol = 0.001, component 1)
qqnorm(resid(modelP))

# overdispersed? negative binomial more appropriate

# best model is modelP, including morphotype, SLA and P, with species specific response to P and site also as random effect
modelSLA=glmer.nb(value ~meanSLA.standard+percP.standard+(1+meanSLA.standard|variable)+(1|site),  data    =hem.both.long)
modelP=glmer.nb(value ~meanSLA.standard+percP.standard+(1+percP.standard|variable)+(1|site),  data    =hem.both.long)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0279573
## (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0294294
## (tol = 0.001, component 1)
modelN=glmer.nb(value ~percN.standard+percP.standard+(1+percN.standard|variable)+(1|site),  data    =hem.both.long)

anova(modelP,modelSLA,modelN)
## Data: hem.both.long
## Models:
## modelP: value ~ meanSLA.standard + percP.standard + (1 + percP.standard | 
## modelP:     variable) + (1 | site)
## modelSLA: value ~ meanSLA.standard + percP.standard + (1 + meanSLA.standard | 
## modelSLA:     variable) + (1 | site)
## modelN: value ~ percN.standard + percP.standard + (1 + percN.standard | 
## modelN:     variable) + (1 | site)
##          Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## modelP    8 3009.1 3054.5 -1496.6   2993.1                        
## modelSLA  8 3011.2 3056.6 -1497.6   2995.2     0      0          1
## modelN    8 3015.2 3060.5 -1499.6   2999.2     0      0          1
modelPalone=glmer.nb(value ~percP.standard+(1+percP.standard|variable)+(1|site),  data    =hem.both.long)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0289687
## (tol = 0.001, component 1)
modelPmorph=glmer.nb(value ~morphotype.x+percP.standard+(1+percP.standard|variable)+(1|site),  data    =hem.both.long)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0193233
## (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge: degenerate Hessian with 1
## negative eigenvalues
modelPSLA=glmer.nb(value ~meanSLA.standard+(1+percP.standard|variable)+(1|site),  data    =hem.both.long)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge: degenerate Hessian with 1
## negative eigenvalues
anova(modelP,modelPalone,modelPmorph)
## Data: hem.both.long
## Models:
## modelPalone: value ~ percP.standard + (1 + percP.standard | variable) + (1 | 
## modelPalone:     site)
## modelP: value ~ meanSLA.standard + percP.standard + (1 + percP.standard | 
## modelP:     variable) + (1 | site)
## modelPmorph: value ~ morphotype.x + percP.standard + (1 + percP.standard | 
## modelPmorph:     variable) + (1 | site)
##             Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)
## modelPalone  7 3011.3 3051.0 -1498.7   2997.3                          
## modelP       8 3009.1 3054.5 -1496.6   2993.1  4.1927      1  0.0405987
## modelPmorph  9 2999.5 3050.5 -1490.7   2981.5 11.6537      1  0.0006407
##                
## modelPalone    
## modelP      *  
## modelPmorph ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelwP=glmer.nb(value ~morphotype.x+percP.standard+(1|variable)+(1|site),  data    =hem.both.long)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0253589
## (tol = 0.001, component 1)
modelPwallP=glmer.nb(value ~morphotype.x+(1|variable)+(1|site),  data    =hem.both.long)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0224301
## (tol = 0.001, component 1)
anova(modelPmorph,modelwP,modelPwallP)
## Data: hem.both.long
## Models:
## modelPwallP: value ~ morphotype.x + (1 | variable) + (1 | site)
## modelwP: value ~ morphotype.x + percP.standard + (1 | variable) + (1 | 
## modelwP:     site)
## modelPmorph: value ~ morphotype.x + percP.standard + (1 + percP.standard | 
## modelPmorph:     variable) + (1 | site)
##             Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)
## modelPwallP  6 3009.2 3043.2 -1498.6   2997.2                          
## modelwP      7 3011.1 3050.8 -1498.6   2997.1  0.0485      1  0.8257047
## modelPmorph  9 2999.5 3050.5 -1490.7   2981.5 15.6670      2  0.0003962
##                
## modelPwallP    
## modelwP        
## modelPmorph ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model with morphotype and perc P, and random effect of percP is best model
summary(modelPmorph)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.2625)  ( log )
## Formula: value ~ morphotype.x + percP.standard + (1 + percP.standard |  
##     variable) + (1 | site)
##    Data: hem.both.long
## 
##      AIC      BIC   logLik deviance df.resid 
##   2999.5   3050.5  -1490.7   2981.5     2133 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.5098 -0.3360 -0.2173 -0.1338 11.8310 
## 
## Random effects:
##  Groups   Name           Variance Std.Dev. Corr
##  variable (Intercept)    3.5283   1.8784       
##           percP.standard 0.2258   0.4752   0.14
##  site     (Intercept)    0.1920   0.4382       
## Number of obs: 2142, groups:  variable, 42; site, 4
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -2.51629    0.40967  -6.142 8.13e-10 ***
## morphotype.xG   0.88645    0.22141   4.004 6.24e-05 ***
## morphotype.xP   0.10121    0.28266   0.358     0.72    
## percP.standard -0.01331    0.17780  -0.075     0.94    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mrph.G mrph.P
## morphtyp.xG -0.305              
## morphtyp.xP -0.235  0.126       
## prcP.stndrd  0.043 -0.394  0.337
dotplot(ranef(modelP, condVar=T))
## $variable

## 
## $site

# check for overdispersion:
pcsq<-sum(resid(modelP,type="pearson")^2)
rdf<-df.residual(modelP)
pchisq(pcsq,rdf,lower.tail=F)
## [1] 1
pcsq/rdf
## [1] 0.8416044
library(blmeco)
## Loading required package: MASS
dispersion_glmer(modelP)
## [1] 0.6364843
# alternatively could use package DHArma



################# Figures 2015 and 2014 data

cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000", "green")
ggplot(data=hem.both.long, aes(x=percP.standard, y=log(value), group=variable, colour=variable)) +
  #    geom_point()+
  geom_smooth(method=lm, se=F,level=0.95)+
  ylab("Abundance")+ xlab("standardized perc P")+
    theme_bw()
## Warning: Removed 1795 rows containing non-finite values (stat_smooth).

hemtop10<-hem.both.long[hem.both.long$variable==c("GRE_PSI"),]
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("OCE_VUL"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("OPU_SP"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("LEI_SP"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("ORT_MET"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("KOA_HAW"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("SAR_ADO"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("PAR_PYR"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("TYM_TYM"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("OLI_SP"),])


cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000", "green", "#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000", "green")
formatlines<-c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6)

ggplot(data=hemtop10, aes(x=percP.standard, y=log(value), group=variable, linetype=variable, colour=variable)) +
  geom_smooth(method=lm, se=F,level=0.95)+
  ylab("Abundance")+ xlab("standardized foliar phosphorus %")+
  scale_linetype_manual(values=formatlines, name="Species")+
  scale_colour_manual(values=cbPalette, name="Species")+
  theme_bw()
## Warning: Removed 297 rows containing non-finite values (stat_smooth).

ggplot(data=hemtop10, aes(x=meanSLA.standard, y=log(value), group=variable, linetype=variable, colour=variable)) +
  geom_smooth(method=lm, se=F,level=0.95)+
  ylab("Abundance")+ xlab("standardized SLA%")+
  scale_linetype_manual(values=formatlines, name="Species")+
  scale_colour_manual(values=cbPalette, name="Species")+
  theme_bw()
## Warning: Removed 297 rows containing non-finite values (stat_smooth).

ggplot(data=hemtop10, aes(x=morphotype.x, y=log(value), group=variable, linetype=variable, colour=variable)) +
  geom_point()+
  ylab("Abundance")+ xlab("standardized foliar phosphorus %")+
  scale_shape_manual(values=formatlines, name="Species")+
  scale_colour_manual(values=cbPalette, name="Species")+
  theme_bw()

Conclusions



(i) How do herbivore communities respond to phenotypic variation?

There is variation in species composition across sites and metrosideros traits. NMDS plots indicate differences between sites and volcanoes in species composition.

(ii) Which plant traits account for herbivore community responses to host-plant phenotype?

CCA results point to the importance of foliar nitrogen, leaf morphotype and specific leaf area. However, there is still a large part of variation not explained.

(iii) Which insect traits account for herbivore community responses to host-plant phenotype

Community composition varies across metrosideros traits. Feeding guild is correlated with foliar %N, where leaf chewers (sig) and sap feeders (NS) are more abundant on higher nutrient trees. Feeding guild is correlated with morphotype, where high trichome density morphotype is correlated with higher abundance gallers and seed feeders, and lower abundances of chewers and sap feeders. Morphotype is positively correlated with proportion specialists (pubescent leaves show relatively higher abundances of specialists). Location of the nymph varies little across morphotypes,